
############################################
############################################
# previous: 4-VAP_VarSelect ###############
############################################

###### Summary ############################################################################################
#This script creates .png images showing maps of the current habiat suitabilty predicted by maxent for each modelling unit. Results are shown 
#globally, within 10 decimal degrees of occurrence records, and for several different thresholding options suggested by Maxent. Maps also show 
#current WHO distribution estimates, and which records were or were not used for models based on their location accuracy.
###########################################################################################################

############################################
# next: 6-VAP_FutProj ######################
############################################

#######################################################################################################
##### MaxEnt3 Plotting ################################################################################
#######################################################################################################

##########
#Settings#
##########

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(stringr)
library(sf)

#basic model settings
#set project name:
projectName<-"WHOSnakes"

#define variable selection metrics:
crits<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)

#define directories
MyDir<-paste("/home/jc217070/Maxent",sep="")
SWDDir<-paste(MyDir, "/Maxent_SWDs/spp_swds",sep="")
OCCDir<-paste(MyDir, "/Maxent_SWDs/spp_occs",sep="")
OutDirCrossval<-paste(MyDir,"/Outputs2_Crossval",sep="")
OutDirFinal<-paste(MyDir,"/Outputs3_Final",sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")
dir.create(path=file.path(MyDir,"tmp.pbs"),showWarnings = FALSE)
PbsDir <- paste(MyDir, "/tmp.pbs", sep="")
setwd(PbsDir)

#read in lookup tables
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin3.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax

#read in base layers
worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))
mybg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")

#define plot colors
my_col1 = rev(terrain.colors(n = 20, alpha=1))
my_col2 = rev(terrain.colors(n = 20, alpha=0.5)) #alpha=0:1, higher = less transparent
my_col3 = c("white","green")
my_col4 = c("white","grey", "green")
mycols1<-c("grey90","grey85","grey80","grey75","grey70",
           "grey65","grey60","grey55","grey50","grey45",
           "grey40","grey35","grey30")
mycols2<-c("royalblue","cyan","navy","skyblue","gold",
           "red","darkred","darkorange","orangered",
           "deeppink","pink","purple","slateblue1","lightsalmon","burlywood4","coral3",
           "deeppink4","darkorange4","cyan4","purple4")
WHOcol <- rgb(150, 180, 250, max = 255, alpha = 0) #alpha=0 means completely transparent

##############
# main script#
##############

mapspp<-unique(DatSum$ModUnit)
mapspp<-paste("ModUnit_", gsub(" ","_",mapspp), sep="")

#mapsppB<-unique(DatSum$ModUnit[DatSum$Priority==0 | DatSum$Priority==1])
#mapsppB<-unique(DatSum$ModUnit[DatSum$Redo==1])
#mapsppB<-c("Walterinnesia aegyptia","Walterinnesia morgani")
#mapsppB<-unique(DatSum$ModUnit [DatSum$MergeShort %in% c("Bitis arietans","Bitis gabonica group","Bitis nasicornis",
#                                                          "Cerastes cerastes complex","Echis coloratus super-complex",
#                                                          "Dendroaspis polylepis","Dendroaspis viridis group",
#                                                          "Naja-Afronaja complex","Naja-Uraeus complex",
#                                                          "Naja-Boulengerina annulata group","Naja-Boulengerina melanoleuca group",
#                                                          "Pseudohaje spp.","Thelotornis spp.","Dispholidus typus")])
 
mapsppB<-unique(DatSum$ModUnit)
mapsppB<-paste("ModUnit_", gsub(" ","_",mapsppB), sep="")

for (sp in mapspp) { #for each modelling unit...

  ##########################
  if (sp %in% mapsppB){
  ##########################
  
    
    
  n<-which(mapspp==sp) #gets row number for current species (for reference back to full data frame)
  spX<-gsub("ModUnit_","",sp)
  is<-which(DatSum$ModUnit==gsub("_"," ",spX))
    
  for (crit in crits){
    
  if (file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc.gz", sep="")) |
      file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""))){ #safety, skip species with no output
    
    #print out summary of whats going on (keep track of progress: start)
    Sys.sleep(0.1)
    print(paste(Sys.time(),"plotting", sp,sep=" "))
    
    #get species records used for model:
    sprecs<-read.table(paste(OCCDir,"/",spX,"_mod.csv",sep=""), sep=",", header=TRUE)   #read in species recs
    recnum<-length(sprecs[,1])
    
    # read in the maxent results thresholds
    results <- read.csv(paste(OutDirFinal, "/", crit, "/", sp, "/maxentResults.csv", sep=""))
    TrainAUC<-results$Training.AUC[1]
    resultsCV <- read.csv(paste(OutDirCrossval, "/", crit, "/", sp, "/maxentResults.csv", sep=""))
    finrow<-length(resultsCV[,1])
    TestAUC<-paste(resultsCV$Test.AUC[finrow])
    
    Oms<-read.csv(paste(OutDirFinal, "/", crit, "/", sp, "/",gsub("ModUnit_","",sp),"_omission.csv", sep=""))
    
    Omsn<-which(Oms$Training.omission==max(Oms$Training.omission[Oms$Training.omission<=0.05]))
    myt<-Oms$Corresponding.Cloglog.value[Omsn]
    
    # thresholds
    thresh1<-as.numeric(results[1,c("Maximum.training.sensitivity.plus.specificity.Cloglog.threshold")])
    thresh2<-as.numeric(results[1,c("Equate.entropy.of.thresholded.and.original.distributions.Cloglog.threshold")])
    thresh3<-as.numeric(results[1,c("Balance.training.omission..predicted.area.and.threshold.value.Cloglog.threshold")])
    thresh4<-myt
    
    threshs<-c(thresh1,thresh2,thresh3,thresh4)
    #t<-as.numeric(DatSum$thresh[is])[1] #only use first one as they are all the same for each species within this mapping unit
    t <- 3
    thresh<-threshs[t]

    DatSum$thresh[is]<-thresh #write out threshold used for each species into summary table
    write.table(x=DatSum, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin3.csv",sep=""),   
                na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
    
    if(length(sprecs[,1])>0){
      
      Sys.sleep(0.1)
      print(paste(Sys.time(), "making extent for modelling unit", sp, sep=" "))
      
    
    rangewidthDG<-max((max(na.omit(sprecs$long))-min(na.omit(sprecs$long))),
                      (max(na.omit(sprecs$lat))-min(na.omit(sprecs$lat))))
    rangewidthM<-rangewidthDG*100000 #turn dec deg into meters
    rangewidthM<-ifelse(rangewidthM>3000000,3000000,rangewidthM) #make buffer no bigger than this
    rangewidthM<-ifelse(rangewidthM<1000000,1000000,rangewidthM) #make buffer no smaller than this
    y<-rangewidthM
    
    #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
    xrange<- (max(na.omit(sprecs$long)))-(min(na.omit(sprecs$long)))
    yrange<- (max(na.omit(sprecs$lat)))-(min(na.omit(sprecs$lat)))
    xcenter<- (min(na.omit(sprecs$long)))+(xrange/2)
    ycenter<- (min(na.omit(sprecs$lat)))+(yrange/2)
    rangemax<-max(xrange,yrange)
    xmin<-ifelse((xcenter-(rangemax/2)- y/100000) >= -180, (xcenter-(rangemax/2)- y/100000), -180)
    xmax<-ifelse((xcenter+(rangemax/2)+ y/100000) <=  180, (xcenter+(rangemax/2)+ y/100000),  180)
    ymin<-ifelse((ycenter-(rangemax/2)- y/100000) >=  -90, (ycenter-(rangemax/2)- y/100000),  -90)
    ymax<-ifelse((ycenter+(rangemax/2)+ y/100000) <=   90, (ycenter+(rangemax/2)+ y/100000),   90)
    myext<-extent(xmin,xmax,ymin, ymax)
  } else{
    myext<-extent(-180,180,-90,90)
  }
    
    
    ####################################
    #Layer processing:
    
    Sys.sleep(0.1)
    print(paste(Sys.time(),"making & writing thresholded SDM versions",sep=" "))
    
    print(paste(Sys.time(),"unzipping & reading raw model",sep=" "))
    if(!(file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep="")))){
      gunzip(filename=paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc.gz", sep=""), 
             destname=paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""),
             remove=FALSE, overwrite=TRUE, skip=TRUE)
    }
    SDMraw<-raster(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""), crs="+init=epsg:4326")
    
    print(paste(Sys.time(),"cropping & thresholding for vetting plots",sep=" "))
    mybgxs<-crop(mybg, y=myext)
    SDM0 <- crop(SDMraw, y=myext)
    writeRaster(SDM0, filename=paste(OutDirFinal, "/", crit, "/", sp, "/",sp,"_CurrentCrop.asc",sep=""), format="ascii", overwrite=TRUE)
    
    SDM1 <- SDM0
    SDM2 <- SDM0
    SDM3 <- SDM0
    SDM4 <- SDM0
    values(SDM1)<-ifelse(values(SDM1)<threshs[1], 0, values(SDM1))
    values(SDM2)<-ifelse(values(SDM2)<threshs[2], 0, values(SDM2))
    values(SDM3)<-ifelse(values(SDM3)<threshs[3], 0, values(SDM3))
    values(SDM4)<-ifelse(values(SDM4)<threshs[4], 0, values(SDM4))
    writeRaster(SDM3, filename=paste(OutDirFinal, "/", crit, "/", sp, "/",sp,"_CurrentCropThresh.asc",sep=""), format="ascii", overwrite=TRUE)
    SDMstack<-stack(SDM1,SDM2,SDM3,SDM4)
    
    ######################################
    #Plotting:
    
    Sys.sleep(0.1)
    print(paste(Sys.time(),"plotting results",sep=" "))

    png(filename=paste(PlotDir,"/",sp,"_SDM.png",sep=""), units="in", width=10, height=12, res=500) #png smaller file size hat pdf
    par(mfrow = c(3,2), mar = c(1,1,1,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
    #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)

    ############
    #first plot#
    ############
    
    print(paste(Sys.time(),"plot 1",sep=" "))
    plot(SDMraw, col=my_col1, box=FALSE, bty="n", zlim=c(0,1), axes=FALSE)
    plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
    
    for (mapunit in DatSum[is,"gen_sp_subtax"]){ #for each species unit in modelling unit
      k<-which(DatSum$gen_sp_subtax==mapunit)
      yesCountry<-str_split(DatSum[k,"country"], "/")                   #split countries of occurrence into separate elements
      mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
      plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
      
    if(file.exists(paste(OCCDir,"/",mapunit,"_allRecs.csv",sep=""))){
      urecs<-read.table(paste(OCCDir,"/",mapunit,"_allRecs.csv",sep=""), sep=",", header=TRUE)
      points(urecs$lat~urecs$long, col="black", pch=1, cex=1)
    }
    }

    text(-172,-5, labels=expression("O = all records"),pos=4, col="black", cex=0.8)
    text(-172,-15, labels=expression("+ = model records"),pos=4, col="black", cex=0.8)
    text(-172,-25, labels=paste("Training AUC = ",TrainAUC,sep=""),pos=4, col="black", cex=0.8)
    text(-172,-35, labels=paste("Test AUC = ",TestAUC,sep=""),pos=4, col="black", cex=0.8)
    text(-172,-45, labels=paste("Threshold: T",t, " = ", round(thresh,2),sep=""),pos=4, col="black", cex=0.8)
    
    legend(legend=paste("Raw global suitability model"),x="topright", cex=1.5, inset=c(0.02,0.04), bty="n")
    
    #############
    #second plot#
    #############
    
    print(paste(Sys.time(),"plot 2",sep=" "))
    plot(mybgxs, col=mycols1, ylab="latitude", xlab="longitude", box=FALSE, bty="l", cex.axis=0.8, legend=FALSE)
    plot(SDM0, col=my_col2, legend=FALSE, add=TRUE)
    plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)

    for (mapunit in DatSum[is,"gen_sp_subtax"]){
      k<-which(DatSum$gen_sp_subtax==mapunit)
      q<-which(DatSum[is,"gen_sp_subtax"]==mapunit)
      yesCountry<-str_split(DatSum[k,"country"], "/")                   #split countries of occurrence into separate elements
      mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
      plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
    }
    
    for (mapunit in DatSum[is,"gen_sp_subtax"]){
      k<-which(DatSum$gen_sp_subtax==mapunit)
      q<-which(DatSum[is,"gen_sp_subtax"]==mapunit)
      species_shapeA <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",paste(gsub("ModUnit_","",mapunit))))
      species_shapeB <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",paste(DatSum$WHOSpecies[k])))
      
      if(length(species_shapeA)>=1){
        plot(species_shapeA[2], add=TRUE, border=mycols2[q], col=WHOcol, lwd= 1)
      }else{
        if(length(species_shapeB)>=1){
          plot(species_shapeB[2], add=TRUE, border=mycols2[q], col=WHOcol, lwd= 1)
        }
      }
      
      if(file.exists(paste(OCCDir,"/",mapunit,".csv",sep=""))){
        urecs<-read.table(paste(OCCDir,"/",mapunit,"_allRecs.csv",sep=""), sep=",", header=TRUE)
      #  points(urecs$lat~urecs$long, col=mycols2[q], pch=1, cex=1)
      }
      if(file.exists(paste(SWDDir,"/",mapunit,".csv",sep=""))){
        urecs2<-read.table(paste(SWDDir,"/",mapunit,".csv",sep=""), sep=",", header=TRUE)
      #  points(urecs2$lat~urecs2$long, col=mycols2[q], pch=3, cex=1)
      }
      
      if(file.exists(paste(OCCDir,"/",mapunit,".csv",sep=""))){
        if(length(urecs[,1])>=1){
          myx<-mean(urecs$long)
          myy<-min(urecs$lat)-2
          graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[q])
        }
      }
    }
    legend(legend=paste("T0= no threshold"),x="topright", cex=1.5, inset=c(0.02,0.02),bg="white")
    
    ###############
    #plots 3,4,5,6#
    ###############
    
    for (p in c(1:4)){
      print(paste(Sys.time(),"plot",p+2 ,sep=" "))
      plot(mybgxs, col=mycols1, ylab="latitude", xlab="longitude", box=FALSE, bty="l", cex.axis=0.8, legend=FALSE)
      plot(SDMstack[[p]], col=my_col2, legend=FALSE, add=TRUE)
      plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
      for (mapunit in DatSum[is,"gen_sp_subtax"]){
        k<-which(DatSum$gen_sp_subtax==mapunit)
        q<-which(DatSum[is,"gen_sp_subtax"]==mapunit)
        yesCountry<-str_split(DatSum[k,"country"], "/")                   #split countries of occurrence into separate elements
        mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
        plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
        if(file.exists(paste(OCCDir,"/",mapunit,".csv",sep=""))){
          urecs<-read.table(paste(OCCDir,"/",mapunit,"_allRecs.csv",sep=""), sep=",", header=TRUE)
          points(urecs$lat~urecs$long, col=mycols2[q], pch=1, cex=1)
          }
        if(file.exists(paste(SWDDir,"/",mapunit,".csv",sep=""))){
          urecs2<-read.table(paste(SWDDir,"/",mapunit,".csv",sep=""), sep=",", header=TRUE)
          points(urecs2$lat~urecs2$long, col=mycols2[q], pch=3, cex=1)
        }
        }
      legend(legend=paste("T",p,"=", round(threshs[p],2),sep=""),x="topright", cex=1.5, inset=c(0.02,0.02),bg="white")
    }
    
    ############
    #plot title#
    ############
    
    mtext(paste("'",gsub("_"," ",spX),"'"," modelling unit",sep=""), side = 3, line = 1, outer = TRUE, cex= 1.5)
    
    dev.off()
    
    #print out summary of whats going on (keep track of progress: finish)
    Sys.sleep(0.1)
    print(paste(Sys.time(),"done plotting", sp,sep=" "))

    if((file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc.gz", sep="")))){
      file.remove(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""))
    } else {
      gzip(filename=paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""),overwrite=TRUE, remove=TRUE)
    }
    
    #if((file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_limitFactor.asc", sep="")))){
    #  gzip(filename=paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_limitFactor", sep=""),overwrite=TRUE, remove=TRUE)
    #}
    
    gzip(filename=paste(OutDirFinal, "/", crit, "/", sp, "/",sp,"_CurrentCrop.asc",sep=""),overwrite=TRUE, remove=TRUE)
    gzip(filename=paste(OutDirFinal, "/", crit, "/", sp, "/",sp,"_CurrentCropThresh.asc",sep=""),overwrite=TRUE, remove=TRUE)
    
  } #end for safety 
  }
  }
} #end for species loop


############################################
############################################
# next: 6-VAP_FutProj ######################
############################################


